%function [sig_diff_ttest, l_c_mean, l_c_ste, l_e_mean, l_e_ste, cells_total, cells_good, c_mat, e_mat]=analyze_hom_data(hom_data,layer, thresh, integ)




% condition: -2 -1 1 2, neg: dark, 2: running
cond=[-2];
% cortical layer: 2 or 5
layer=[2 5];
% animals 0: all 1:select
animals=[1:9 10 11 12 13];
% cell types: 99:cells or 100:dendrites
types=[99];
% use active cells only
use_act_only=0;
% normalize to first 2 datapoints
norm_to_baseline=0;
% count active cells instead of taking mean activity
count_act_cells=0;

% %integ is 1 for integration method, 0 for cell counting method
integ=1;
%
% %threshold
thresh=1.1;

% using raw data
raw_data_flag=1;


% security idiot with flashlight
hom_data(43).raw_data(6).act_code(7000:11000)=0;


cells_good=[];
cells_total=[];


cnt=0;
cnt_e=[];
cnt_c=[];
act_lp_corr=[];
e_act={[] [] [] [] [] [] []};
c_act={[] [] [] [] [] [] []};

tpoints=[-48 -24 6 18 24 48 72];

%tpoints=[-48 -24 12 24 48 72];

e_mat=cell(size(hom_data,2),7);
c_mat=cell(size(hom_data,2),7);

vars_e={'lp_e','mean_templ_e','stim_dur_e','running_act_e'};
vars_c={'lp_c','mean_templ_c','stim_dur_c','running_act_c'};

for ind=1:length(vars_e)
    eval([vars_e{ind} '=-1*ones(length(hom_data),7);']);
end
for ind=1:length(vars_c)
    eval([vars_c{ind} '=-1*ones(length(hom_data),7);']);
end

all_act=[];
for ind=1:length(hom_data)
    act_cell_ind=ones(1,length(sum(hom_data(ind).raw_data(1).sig_act)));
    for knd=1:length(hom_data(ind).raw_data)
        act_cell_ind=act_cell_ind&sum(hom_data(ind).raw_data(knd).sig_act)>10;
    end
    
    
    for knd=1:length(hom_data(ind).raw_data)
        curr_ind=find(hom_data(ind).raw_data(knd).timepoint==tpoints);
        l_ind=zeros(length(hom_data(ind).raw_data(knd).act_code),1);
        for jnd=1:length(cond)
            l_ind=l_ind+double(hom_data(ind).raw_data(knd).act_code==cond(jnd));
        end
        l_ind=logical(l_ind);
        
        cell_ind=zeros(1,size(hom_data(ind).raw_data(knd).sig_act,2));
        for jnd=1:length(types)
            cell_ind=cell_ind+double(hom_data(ind).raw_data(knd).ROItype==types(jnd));
        end
        %
        if knd==1
            cells_total(ind)=sum(cell_ind);
        end
        cell_ind=logical(cell_ind);
        
        fcells=filled_cells_all;
        if fcells{ind}~=-1
            cell_ind(fcells{ind})=0;
        end
        
        if knd==1
            cells_good(ind)=sum(cell_ind);
        end
        cell_ind=logical(cell_ind);
        
        
        if use_act_only
            cell_ind=cell_ind&act_cell_ind;
        end
        if count_act_cells
            tmp=sum(hom_data(ind).raw_data(knd).sig_act(l_ind,cell_ind))>10/sum(l_ind);
        elseif raw_data_flag
            tmp=hom_data(ind).raw_data(knd).sig_act(l_ind,cell_ind);
            tmp=tmp./repmat(median(tmp),size(tmp,1),1);
            
            
            if integ==1  %integral method then integ=1, cell sum method integ=0
                %
                tmp(tmp<thresh)=0;
                tmp=sum(tmp)/sum(l_ind);
            elseif integ==2
                tmp=max(tmp);
            else
                % % %                 for jnd=1:size(tmp,2)
                % % %                     std_noise=sqrt(sum((tmp(tmp(:,jnd)<1,jnd)-1).^2/sum(tmp(:,jnd)<1)));
                % % %                     tmp(:,jnd)=tmp(:,jnd)>(1+10*std_noise);
                % % %                 end
                % % %                 tmp=sum(tmp)/sum(l_ind);
                
                
                tmp=sum(tmp>thresh)/sum(l_ind);
                all_act(end+1:end+length(tmp))=tmp;
                tmp=tmp>0;
            end
        else
            tmp=sum(hom_data(ind).raw_data(knd).sig_act(l_ind,cell_ind))/sum(l_ind);
        end
        if sum(hom_data(ind).animal_ID==animals)
            if sum(hom_data(ind).layer==layer)
                if hom_data(ind).animal_ID<=9
                    sum(tmp)/length(tmp);
                    cnt_e=[cnt_e ind];
                    e_act{curr_ind}(end+1:end+length(tmp))=tmp;
                    e_mat{ind,curr_ind}(end+1:end+length(tmp))=tmp;
                    mean_templ_e(ind,curr_ind)=mean(hom_data(ind).raw_data(knd).template(:));
                    stim_dur_e(ind,curr_ind)=sum(l_ind);
                    running_act_e(ind,curr_ind)=sum((hom_data(ind).raw_data(knd).act_code)==-2)/sum((hom_data(ind).raw_data(knd).act_code)<0);
                else
                    cnt_c=[cnt_c ind];
                    c_act{curr_ind}(end+1:end+length(tmp))=tmp;
                    c_mat{ind,curr_ind}(end+1:end+length(tmp))=tmp;
                    mean_templ_c(ind,curr_ind)=mean(hom_data(ind).raw_data(knd).template(:));
                    stim_dur_c(ind,curr_ind)=sum(l_ind);
                    running_act_c(ind,curr_ind)=sum((hom_data(ind).raw_data(knd).act_code)==-2)/sum((hom_data(ind).raw_data(knd).act_code)<0);
                end
                cnt=cnt+1;
            end
        end
    end
end


% for ind=1:length(vars_e)
%     eval([vars_e{ind} '=' vars_e{ind} '(unique(cnt_e),:);']);
%     eval([vars_e{ind} '(' vars_e{ind} '==-1)=NaN;']);
%     if norm_to_baseline
%         for knd=1:size(e_mat,1)
%             eval([vars_e{ind} '(knd,:)=' vars_e{ind} '(knd,:)/nanmean(' vars_e{ind} '(knd,1:2));']);
%         end
%     end
% end
%
% for ind=1:length(vars_c)
%     eval([vars_c{ind} '=' vars_c{ind} '(unique(cnt_c),:);']);
%     eval([vars_c{ind} '(' vars_c{ind} '==-1)=NaN;']);
%     if norm_to_baseline
%         for knd=1:size(c_mat,1)
%             eval([vars_c{ind} '(knd,:)=' vars_c{ind} '(knd,:)/nanmean(' vars_c{ind} '(knd,1:2));']);
%         end
%     end
% end

% e_mat=e_mat(unique(cnt_e),:);
% c_mat=c_mat(unique(cnt_c),:);
% lp_e=lp_e(unique(cnt_e),:);
% lp_c=lp_c(unique(cnt_c),:);
% mean_templ_e=mean_templ_e(unique(cnt_e),:);
% mean_templ_c=mean_templ_c(unique(cnt_c),:);
%
% e_mat(e_mat==-1)=NaN;
% c_mat(c_mat==-1)=NaN;
% lp_e(lp_e==-1)=NaN;
% lp_c(lp_c==-1)=NaN;
% mean_templ_e(mean_templ_e==-1)=NaN;
% mean_templ_c(mean_templ_c==-1)=NaN;
%
% for ind=1:size(e_mat,1)
%     e_mat(ind,:)=e_mat(ind,:)/nanmean(e_mat(ind,1:2));
%     lp_e(ind,:)=lp_e(ind,:)/nanmean(lp_e(ind,1:2));
%     mean_templ_e(ind,:)=mean_templ_e(ind,:)/nanmean(mean_templ_e(ind,1:2));
% end
% for ind=1:size(c_mat,1)
%     c_mat(ind,:)=c_mat(ind,:)/nanmean(c_mat(ind,1:2));
%     lp_c(ind,:)=lp_c(ind,:)/nanmean(lp_c(ind,1:2));
%     mean_templ_c(ind,:)=mean_templ_c(ind,:)/nanmean(mean_templ_c(ind,1:2));
% end
l_c_mean=[];
l_e_mean=[];
l_e_ste=[];
l_c_ste=[];

for ind=1:length(e_act)
    if norm_to_baseline
        ee=e_act{ind}/nanmean([nanmean(e_act{1}) nanmean(e_act{2})]);
        cc=c_act{ind}/nanmean([nanmean(c_act{1}) nanmean(c_act{2})]);
    else
        ee=e_act{ind};
        cc=c_act{ind};
    end

    l_e_mean(ind)=nanmean(ee);
    l_c_mean(ind)=nanmean(cc);
    l_e_ste(ind)=nanstd(ee)/sqrt(sum(~isnan(ee)));
    l_c_ste(ind)=nanstd(cc)/sqrt(sum(~isnan(cc)));
    try
        [~,sig_diff_ttest(ind)]=ttest2(ee,cc);
        sig_diff_rsum(ind)=ranksum(ee,cc);
    catch
    end
end


figure;
clf
hold on
rgb_val=[1 1 1]*0;
plot(tpoints(1:6),l_c_mean(1:6)-l_e_mean(1:6),'o','color',rgb_val,'linewidth',2)
plot(tpoints(1:2),l_c_mean(1:2)-l_e_mean(1:2),':','color',rgb_val,'linewidth',2)
plot(tpoints(3:6),l_c_mean(3:6)-l_e_mean(3:6),':','color',rgb_val,'linewidth',2)
for ind=1:6
    plot([1 1]*tpoints(ind),[-1 1]*l_e_ste(ind)+l_c_mean(ind)-l_e_mean(ind),'color',rgb_val)
end

figure;
clf
hold on
rgb_val=[0 0 0]*1;
plot(tpoints(1:6),l_e_mean(1:6),'o','color',rgb_val,'linewidth',2)
plot(tpoints(1:2),l_e_mean(1:2),'color',rgb_val,'linewidth',2)
plot(tpoints(3:6),l_e_mean(3:6),'color',rgb_val,'linewidth',2)
for ind=1:6
    plot([1 1]*tpoints(ind),[-1 1]*l_e_ste(ind)+l_e_mean(ind),'color',rgb_val)
end

plot(tpoints(1:6),l_c_mean(1:6),'or','linewidth',2)
plot(tpoints(1:2),l_c_mean(1:2),'r','linewidth',2)
plot(tpoints(3:6),l_c_mean(3:6),'r','linewidth',2)
for ind=1:6
    plot([1 1]*tpoints(ind),[-1 1]*l_c_ste(ind)+l_c_mean(ind),'r')
end

set(gca,'xtick',[tpoints(1:2) 0 tpoints(3:6)])
box off

title_str=({sprintf('Cond: %s - animals; %s',num2str(cond),num2str(animals)) ...
    sprintf('ne: %d, nc %d, integ: %d, thresh: %d',length(e_act{3}),length(c_act{3}),integ,thresh)});
title(title_str,'fontweight','bold');
xlabel('Mean over cells')

e_sites=[];
c_sites=[];


for ind=1:length(hom_data)
    for knd=1:6
        e_sites(ind,knd)=mean(e_mat{ind,knd});
        c_sites(ind,knd)=mean(c_mat{ind,knd});
    end
end
% else
%      for ind=1:4:length(hom_data)
%         for knd=1:6
%             e_sites((ind-1)/4+1,knd)=mean(cell2mat(e_mat(ind:ind+3,knd)'));
%             c_sites((ind-1)/4+1,knd)=mean(cell2mat(c_mat(ind:ind+3,knd)'));
%         end
%     end
% end

e_sites=e_sites/nanmean(nanmean(e_sites(:,1:2)));
m_e_sites=nanmean(e_sites);
s_e_sites=nanstd(e_sites)/sqrt(size(e_sites,1));
c_sites=c_sites/nanmean(nanmean(c_sites(:,1:2)));
m_c_sites=nanmean(c_sites);
s_c_sites=nanstd(c_sites)/sqrt(size(c_sites,1));


% figure;
% clf
% hold on
% plot(tpoints(1:6),m_e_sites(1:6),'ok','linewidth',2)
% plot(tpoints(1:2),m_e_sites(1:2),'k','linewidth',2)
% plot(tpoints(3:6),m_e_sites(3:6),'k','linewidth',2)
% for ind=1:6
%     plot([1 1]*tpoints(ind),[-1 1]*s_e_sites(ind)+m_e_sites(ind),'k')
% end
%
% plot(tpoints(1:6),m_c_sites(1:6),'or','linewidth',2)
% plot(tpoints(1:2),m_c_sites(1:2),'r','linewidth',2)
% plot(tpoints(3:6),m_c_sites(3:6),'r','linewidth',2)
% for ind=1:6
%     plot([1 1]*tpoints(ind),[-1 1]*s_c_sites(ind)+m_c_sites(ind),'r')
% end
% title(title_str,'fontweight','bold');
% xlabel('Mean over regions')
%
%
% % % for ind=1:length(hom_data)
% % %     for knd=1:length(hom_data(ind).raw_data)
% % %         figure(knd);
% % %         clf;
% % %         imagesc(hom_data(ind).raw_data(knd).sig_act');
% % %         title(sprintf('%d - %d - %d', hom_data(ind).animal_ID,hom_data(ind).region_ID,knd ));
% % %     end
% % %     pause
% % % end
%
%%% plot fract running
% fract_running=nan(length(hom_data),7);
% for ind=1:length(hom_data)
%     for knd=1:length(hom_data(ind).raw_data)
%         curr_ind=find(tpoints==hom_data(ind).raw_data(knd).timepoint);
%         fract_running(ind,curr_ind)=sum((hom_data(ind).raw_data(knd).act_code)==2)/sum((hom_data(ind).raw_data(knd).act_code)>0);
%     end
% end

% running_act_c(running_act_c==-1)=nan;
% running_act_e(running_act_e==-1)=nan;
% 
% 
% figure;
% hold on
% plot(nanmean(running_act_e),'k')
% plot(nanmean(running_act_c),'r')









